In [96]:
from __future__ import division, print_function
import os
# Third-party
import matplotlib.pyplot as plt
import numpy as np
%matplotlib inline
# Custom
import gary.dynamics as gd
import gary.integrate as gi
import gary.io as io
import gary.potential as gp
from gary.units import galactic
from streammorphology import potential_registry, read_allfreqs
In [97]:
potential = potential_registry['via-lactea-log']
Read in initial conditions
In [40]:
w0 = np.loadtxt("/Users/adrian/projects/morphology/misc/progenitor_ICs.txt")
w0[:,3:6] = (w0[:,3:6]*u.km/u.s).decompose(galactic).value
classification = w0[:,6]
w0 = w0[:,:6]
Make a directory for freqmapping and save initial conditions with correct units
In [5]:
path = os.path.join("/Users/adrian/projects/morphology/output/ViaLacteaStreams/")
if not os.path.exists(path):
os.mkdir(path)
In [8]:
w0_filename = os.path.join(path, "w0.npy")
np.save(w0_filename, w0)
In [41]:
t,all_orbits = potential.integrate_orbit(w0, dt=2., nsteps=25000,
Integrator=gi.DOPRI853Integrator)
In [42]:
per = np.sqrt(np.sum(all_orbits[:,:,:3]**2,axis=-1)).min(axis=0)
apo = np.sqrt(np.sum(all_orbits[:,:,:3]**2,axis=-1)).max(axis=0)
In [92]:
t1,w1 = potential.integrate_orbit(w0, dt=0.5, nsteps=100000,
Integrator=gi.DOPRI853Integrator)
In [93]:
naff = gd.NAFF(t1[:len(t1)//2+1])
In [94]:
i = 0
fs = [(w1[:len(t1)//2+1,i,j] + 1j*w1[:len(t1)//2+1,i,j+3]) for j in range(3)]
fxyz1,d,ixes = naff.find_fundamental_frequencies(fs, nintvec=15)
fs = [(w1[len(t1)//2:,i,j] + 1j*w1[len(t1)//2:,i,j+3]) for j in range(3)]
fxyz2,d,ixes = naff.find_fundamental_frequencies(fs, nintvec=15)
In [95]:
(fxyz2 - fxyz1) / fxyz1n
Out[95]:
In [98]:
d = read_allfreqs(os.path.join(path, 'allfreqs.dat'), len(w0))
d.dtype.names
Out[98]:
In [100]:
loop = d[d['loop']]
boxy = d[~d['loop']]
In [111]:
loop_nan = np.any(np.any(np.isnan(loop['fRphiz']), axis=-1), axis=-1)
boxy_nan = np.any(np.any(np.isnan(boxy['fxyz']), axis=-1), axis=-1)
In [116]:
loop[loop_nan]['dt'], loop[loop_nan]['nsteps']
Out[116]:
In [117]:
boxy[boxy_nan]['dt'], boxy[boxy_nan]['nsteps']
Out[117]:
Hm, ok, from the above, it looks like these orbits weren't integrated for long enough?
Plot the frequency diffusion rate vs. Emily's by-eye classification, colored by apocenter distance
In [118]:
boxy_fdiff = (boxy['fxyz'][:,1] - boxy['fxyz'][:,0]) / boxy['fxyz'][:,0] / (boxy['dt']*boxy['nsteps']/2.)[:,np.newaxis]
boxy_fdiff = np.max(boxy_fdiff, axis=1)
# boxy_fdiff = boxy_fdiff[np.isfinite(boxy_fdiff)]
plt.figure(figsize=(10,8))
plt.scatter(classification[~d['loop']], boxy_fdiff, c=apo[~d['loop']])
plt.colorbar()
plt.yscale('log')
plt.xlim(-1.,3.)
plt.ylim(1E-14, 1E-3)
Out[118]:
Ok, clearly I should only take apo < 50 kpc or something
In [145]:
ix = apo[~d['loop']] < 45
boxy_fdiff = (boxy['fxyz'][ix,1] - boxy['fxyz'][ix,0]) / boxy['fxyz'][ix,0] / (boxy['dt']*boxy['nsteps']/2.)[ix,np.newaxis]
boxy_fdiff = np.max(boxy_fdiff, axis=1)
loop_fdiff = (loop['fRphiz'][ix,1] - loop['fRphiz'][ix,0]) / loop['fRphiz'][ix,0] / (loop['dt']*loop['nsteps']/2.)[ix,np.newaxis]
loop_fdiff = np.max(loop_fdiff, axis=1)
# for coloring points
boxy_maxfreq = np.max(np.abs(boxy['fxyz'][ix,0]),axis=-1)
loop_maxfreq = np.max(np.abs(loop['fRphiz'][ix,0]),axis=-1)
fig,axes = plt.subplots(1,2,figsize=(12,6),sharex=True,sharey=True)
ax = axes[0]
c = ax.scatter(classification[~d['loop']][ix], boxy_fdiff,
marker='s', alpha=0.5, s=35, c='k')# c=per[~d['loop']][ix]) # c=boxy_maxfreq)
# fig.colorbar(c)
ax.set_yscale('log')
ax.set_xlim(-0.5,2.5)
ax.set_ylim(1E-14, 1E-3)
ax.xaxis.set_ticks([0,1,2])
ax.xaxis.set_ticklabels(['Normal','Fanned','Crazy'])
ax.set_title("Boxy", fontsize=20)
ax = axes[1]
c = ax.scatter(classification[d['loop']][ix], loop_fdiff,
marker='s', alpha=0.5, s=35, c='k') # c=per[d['loop']][ix]) # c=loop_maxfreq)
# fig.colorbar(c)
ax.set_yscale('log')
ax.set_xlim(-0.5,2.5)
ax.set_ylim(1E-14, 1E-3)
ax.xaxis.set_ticks([0,1,2])
ax.xaxis.set_ticklabels(['Normal','Fanned','Crazy'])
ax.set_title("Loop", fontsize=20)
Out[145]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [48]:
loop_fdiff = (loop['fRphiz'][:,1] - loop['fRphiz'][:,0]) / loop['fRphiz'][:,0] / (loop['dt']*loop['nsteps']/2.)[:,np.newaxis]
loop_fdiff = np.max(loop_fdiff, axis=1)
# loop_fdiff = loop_fdiff[np.isfinite(loop_fdiff)]
plt.figure(figsize=(10,8))
plt.scatter(classification[d['loop']], loop_fdiff, c=apo[d['loop']])
plt.colorbar()
plt.yscale('log')
plt.xlim(-1.,3.)
plt.ylim(1E-14, 1E-3)
plt.figure(figsize=(10,8))
plt.scatter(classification[d['loop']], loop_fdiff, c=per[d['loop']])
plt.colorbar()
plt.yscale('log')
plt.xlim(-1.,3.)
plt.ylim(1E-14, 1E-3)
Out[48]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [66]:
these_w0 = w0[w0[:,6] == 1.,:6][:10]
In [67]:
nsteps = 120000
t,w = potential.integrate_orbit(these_w0, dt=1., nsteps=nsteps,
Integrator=gi.DOPRI853Integrator)
In [68]:
for i in range(w.shape[1]):
fig = gd.plot_orbits(w[:nsteps//2], ix=i, marker='.', linestyle='none', alpha=0.1)
fig = gd.plot_orbits(w[nsteps//2:], ix=i, marker='.', linestyle='none', alpha=0.1, axes=fig.axes)
In [69]:
for i in range(w.shape[1]):
E = potential.total_energy(w[:,i,:3], w[:,i,3:6])
dE = np.abs(E[1:] - E[0])
print(dE.max())
In [62]:
naff = gd.NAFF(t[:len(t)//2+1])
fdiff = []
for i in range(w.shape[1]):
fs = [(w[:len(t)//2+1,i,j] + 1j*w[:len(t)//2+1,i,j+3]) for j in range(3)]
fxyz1,d,ixes = naff.find_fundamental_frequencies(fs, nintvec=15)
if (t.max() / np.abs(2*np.pi/fxyz1).max()) < 100:
print("Skipping {}".format(i))
continue
fs = [(w[len(t)//2:,i,j] + 1j*w[len(t)//2:,i,j+3]) for j in range(3)]
fxyz2,d,ixes = naff.find_fundamental_frequencies(fs, nintvec=15)
# diffusion per orbit
fdiff.append(np.abs((fxyz2 - fxyz1) / fxyz1) / (t.max() / np.abs(2*np.pi/fxyz1).max() / 2.))
In [65]:
np.array(map(list, fdiff)).max(axis=1), np.array(map(list, fdiff)).max(axis=1).mean()
Out[65]:
Fractional change of 1E-3 per orbit for 'crazy'
In [63]:
t.max() / np.abs(2*np.pi/fxyz1).max() / 2.
Out[63]:
In [ ]: